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We consider the double exchange model at very low densities. The conditions for the formation 
of self-trapped magnetic polarons are analyzed using an independent polaron model. The issue of 
phase separation in the low density region of the temperature-density phase diagram is discussed. 
We show how electrostatic and localization effects can lead to the substantial suppression of the 
phase separated regime. By examining connections between the resulting phase and the polaronic 
phase, we conclude that they reflect essentially the same physical situation of a ferromagnetic droplet 
containing one single electron. In the ultra diluted regime, we explore the possible stabilization of a 
Wigner crystal of magnetic polarons. Our results are compared with the experimental evidence for 
a polaronic phase in europium hexaboride (EuBg), and we are able to reproduce the experimental 
region of stability of the polaronic phase. We further demonstrate that phase-separation is a general 
feature expected in metallic ferromagnets whose bandwidth depends on the magnetization. 
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I. INTRODUCTION 

The earliest descriptions of the concept of magnetic 
polaron appear with the considerations by de Gennes on 
the relevance of the double exchange (DE) mechanism to 
the mixed valent manganites^, and then with Nagaev's 
studies of antiferromagnetic semiconductors, who coined 
the term ferron also often used in the context^'^. Some 
of the first theories based upon the presence of magnetic 
polarons were developed later in the context of the fer- 
romagnetic semiconductor EuO, to explain the spectacu- 
lar metal-insulator transition found at the onset of ferro- 
magnetism for the Eu-rich samples of this compound'*'^. 
Further developments on this concept permitted the ex- 
planation of the spin-flip Raman scattering characteristic 
of certain diluted magnetic semiconductors (DMS)^'^'^'* 
and, more recently, their presence was also claimed to be 
present in the EuBg hexaborides^'^°. 

Just as in the analogous case of the electrostatic po- 
laron, one can devise a pictorial description of a magnetic 
polaron in real space, consisting of a charge carrier sur- 
rounded by a cloud of polarized local spins in an unpo- 
larized magnetic background — a state that arises from 
the exchange interaction between the carrier spin and the 
lattice spins. A distinction is usually made between the 
so-called bound magnetic polaron (BMP) and the free 
magnetic polaron (FMP). The BMP is invoked when the 
charge carrier is bound via Coulomb interaction to an 
impurity center, and is typical in the magnetic semicon- 
ductors. In this case the trapped carrier polarizes the 
lattice spins within its effective Bohr radius as a conse- 
quence of the s-rf-like interaction between electron and 



lattice spins. The FMP, by contrast, results from the 
fact that a free carrier interacting with lattice spins via 
an s-d coupling, can minimize its kinetic energy by po- 
larizing its vicinity. Under certain conditions this carrier 
can then become self-trapped in the resulting potential 
well created by the effect of the local ferromagnetism. 

The presence of such type of entities is believed to play 
an important role in the emergence of many interest- 
ing properties of several important magnetic materials: 
many peculiarities of the manganites, such as the colos- 
sal magnetoresistance (CMR) effect and other anomalies 
in its transport and magnetism, have been attributed in 
part to the development of magnetic polarons near the 
ferromagnetic transition (although in this case the influ- 
ence of orbital and lattice degrees of freedom is equally 
important)-'^-'^'-'^^'-'^^; an interpretation for the CMR in the 
Mn pyrochlores was also proposed on the basis of mag- 
netic polarons^^; besides their relevance for the already 
mentioned Eu-chalcogenides and II- VI semiconductors, 
in the III-V DMS like Gai-a^Mn^^As, the ferromagnetic 
transition can be interpreted within a BMP percolat- 
ing scenario^^; in EuBg, our target magnetic metal ex- 
hibiting CMR, their presence in signaled in the optical 
rcsponse^'^°, although the theoretical interpretation of 
these results has been subject to questioning^^. 

Given that these different classes of magnetic ma- 
terials have attracted considerable attention in recent 
years because of their potential for the development of 
new magnetoelectronic devices, and since magnetic po- 
larons have an apparently ubiquitous presence among 
them, a number of theoretical approaches to the prob- 
lem have been developed through the years. In particu- 
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lar, extensive work has been done with emphasis in the 
physics of magnetic semiconductors"^'^'^'^^'^* and CMR 
manganites^^'^'''^^'^^. 

In the ensuing sections we will focus our attention on 
the stability conditions for the free magnetic polaron in 
the double exchange model (DEM), with an eye on the 
experimental evidence for magnetic polarons in EuBe^'^*^, 
and having in mind the description of the magneto- 
electronic properties of this material in terms of DE^'^. 
Studies devoted to the polaronic stability in this particu- 
lar model and its variations have been performed by sev- 
eral authors both analytically and numerically, although 
under different assumptions^^'^^'^^'^^'^'''^®'^^'^°. We show 
that the DE-based interpretation of magnetotransport 
in EuBg is consistent with the experimental evidence for 
a polaronic phase mediating the PM-FM transition. In 
particular, within an independent polaron model (IPM) 
we reproduce the experimental temperature and density 
range of the polaronic phases, without adjusting param- 
eters. 

Given that we are focusing on the low density regime 
of the DEM, another issue becomes pressing: the known 
tendency for this model to exhibit a phase separation 
instability at reduced densities'^^''^^'^^''^''''^^''^^. We char- 
acterize the phase-separated regime and study how the 
introduction of electrostatic corrections leads to consider- 
able shrinking of the phase separation region in the phase 
diagram. We discuss how the phase-separated regime is 
connected and compatible with the polaronic phase. 

This paper is organized in the following manner. In 
Sec. II we briefly introduce some experimental details 
regarding electronic transport and magnetism in EuBg, 
only to the extent of motivating our studies. In Sec. IV 
we introduce and discuss an IPM for the DEM at very 
low densities, whose results are then confronted in Sec. IV 
with the experimental evidence for magnetic polarons in 
EuBg. We then address the problem of phase separa- 
tion in the low-density DEM in Sec. V: we discuss the 
emergence of phase separation at low densities, study its 
suppression when electrostatic corrections are taken into 
consideration and explore the connections of the resulting 
phase with the polaronic phase. We consider the ultra 
diluted regime in Sec. VI, and provide estimates for the 
stability of a Wigner crystal of magnetic polarons. In 
Sec. VII we close this paper with our conclusions. 



II. OVERVIEW OF EUBe PROPERTIES 

We briefly review the physics of EuBg that motivates 
our choice of microscopic model. An extensive review of 
the phenomenology of this material was given in Ref. 39, 
and here we limit ourselves to the essential aspects rel- 
evant in the context of this paper. EuBg is a magnetic 
metal^°'^^, with residual resistivities at T — > of the or- 
der of ~ 10/ii7.cm, or less, and exibiths a clear CMR sig- 
nal close to Tc"^'''^ [Fig. 1(a)]. Notwithstanding, EuBg 
has a very small carrier density. More precisely, car- 
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FIG. 1: Some experimental signatures of EuBe. Panel (a) 
shows the behavior of the DC resistivity and magnetic suscep- 
tibility as a function of temperature, as reproduced in Ref. 37. 
Notice the upturn followed by a massive drop in p{T) at Tc, 
which correlates with the onset of magnetic order. In (b) we 
reproduce the results from Ref. 38 for the optical reflectiv- 
ity at different temperatures. The enhancement of ujp as the 
system becomes magnetic is clear. 



rier densities estimated from Hall effect measurements 
amount to as little as ~ lO^^cm"'^, or ~ 0.001 carriers per 
unit celP^''*°'^^. These conduction electrons are believed 
to arise from the presence of defects in the structural 
arrangement of the boron framework*^'^^. Such defects 
generate a surplus of electrons that occupy states in the 
conduction band. 

One of the most intriguing and peculiar features of 
EuBg is arguably the giant and rather unusual blue-shift 
of the unscreened plasma edge, Wp, induced simply by 
a temperature variation^^'^'^ [Fig. 1(b)]. At T > Tc, the 
reflectivity spectrum displays a typical metallic behavior, 
with a very well defined plasma threshold. With the es- 
tablishment of the long-range magnetic order, the plasma 
edge increases markedly in such a way that Wp varies by 
a factor of almost 3 between Tq and T < Tc^^'^^ This 
variation of Wp is consistent with the remarkable enhance- 
ment in the carrier densities as the temperature is lowered 
past Tc^\ 

The conduction and valence bands of EuBg are sepa- 
rated by a gap of the order of 1 eV^^'^^'^-'^'^^. This fun- 
damental gap lies at the X point in the cubic Brillouin 
zone (BZ), and the close proximity of E-pio the bottom 
of the conduction band dictates a pocket-like structure 
for the Fermi surface. Given that the interest will be 
almost completely in a temperature range T < Tq, the 
electronic states in the valence band are disregarded in 
the following. 

The magnetism of EuBg arises entirely from the half- 
filled 4/ shell of Eu2+ in the state ^57/2. This implies 
localized magnetism stemming from magnetic moments 
of magnitude S — 7/2. Within this formulation these 
electrons do not itinerate at all. We designate the re- 
sulting magnetic moment by local spin, and use the term 
magnetization of the system when alluding to long range 
ordered phases of these spins. In addition, for our pur- 
poses the conduction band electrons interact with the lo- 
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cal spins only through the Hund's couphng between the 
electron's spin and the local moments'. 



III. THE INDEPENDENT POLARON MODEL 

The Hamiltonian describing conduction electrons hop- 
ping in a tridimensional cubic lattice, and coupled to 
local spins at each lattice site is 

+ Si ■ fa,(3 cl^Ci^i3 , (1) 

usually known as s — / or Kondo lattice model (KLM) 
Hamiltonian^^ , and is one of the canonical models in cor- 
related electronic systems. In this expression, tij = t 
is the hopping integral between neighboring lattice sites, 
cl ^{ci^a) are the second-quantized fermionic creation (an- 
nihilation) operators at latice site i, Si represents the lo- 
cal magnetic moment of magnitude S = 7/2, Jh the ex- 
change coupling of the latter to the itinerating electrons, 
and f = (ti, r2, ra) is the vector of Pauli matrices. The 
sum in the first term is over all pairs of nearest neighbor- 
ing sites 

Since S ~ 7/2 is a rather high spin, the lo- 
cal spin operator is replaced by the classical vec- 
tor Si , parametrized with spherical angles as Si = 
S(sm{6i)cos{(pi),sm{9i)sm{(pi),cos{9i)y This trans- 
forms the second term in (1) into a generalized poten- 
tial term for the electrons, which will be a disordered 
potential above Tq, where all 5^ are uncorrelated. Fur- 
thermore, the magnetic and electronic time scales are 
well apart in such a way that the magnetic background 
provided by the Si'is essentially quenched. This means 
that the typical time between spin fluctuations is much 
longer than the time for the electronic subsystem to reach 
its ground state 

Under these circumstances, we consider the effective 
DE Hamiltonian^''"''* that obtains in the limit Jh — > oo. 
This is accomplished through a local rotation of the quan- 
tization axis so that it coincides with the direction of 
Si at each site, and projecting out the anti-parallel elec- 
tron states^^. Such anti-parallel states lie higher in en- 
ergy (by ~ Jh) and hence are suppressed in the effective 
Hilbert space. The result is 



Ti-DE 



H. 



(2) 



where the new operators di correspond to an effective 
spinless electron that maintains its spin aligned with each 
local moment, and all information about the magnetic 
background is condensed in the effective hopping ampli- 




FIG. 2: Contrasting the variational approach considered 
in eq. (6) (a) to the more realistic, self-consistent, situation 
(b). Represented are the wave function of the self-localized 
electron, a possible spin configuration and the magnetization 
profile. 
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Within the simple DEM, where no other polaron- 
favoring interactions are included (e.g. not including 
antiferromagnetic (AFM) exchange terms between the lo- 
cal moments.), the only possibility for the stabilization 
of magnetic polarons is at low electronic densities. This 
happens because the local spins of several neighboring 
unit cells are expected to participate in the magnetiza- 
tion cloud of each electron. Were it otherwise (i.e. at high 
electronic densities) , the electronic wave functions would 
overlap considerably destroying this polaron picture. 

In order to investigate this problem we will first dis- 
cuss the thermodynamic stability of magnetic polarons 
within the DEM. We can write a free energy for the 
system including an electronic contribution consisting of 
the ground state energy for the electrons in a given mag- 
netic configuration, the local spins contributing only with 
an entropic term. To tackle the first part, we calculate 
numerically the exact electronic density of states (DOS) 
for a given spin configuration"^^, which we designate by 
{Si}), and extract the disorder-averaged DOS at 
constant magnetization: p{E, M) (M being the normal- 
ized local magnetization) (details about this numerical 
calculation will be given below in Sec. VA). 

With this averaging over disorder, we can write the 
electronic contribution to the free energy as 



Eel(M, Ue) 



J e(^EF{M,ne)-e^Ep{E,M)dE, (4) 



where the dependence of the Fermi energy on both mag- 
netization and electron density was made explicit. When 
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TT-e <C 1; and for the purposes of the current calculatfon, 
Eq. (4) can be approximated simply by 

E,i{M,n,)^ Eb{M)n,, (5) 

with Eb representing the bottom of the band - it becomes 
just a single electron problem. Indeed, given the nature 
of the calculation and approximations involved here, the 
consideration of the finite band filling introduces only 
minor corrections and thus we proceed with the above 
approximation to the electronic energy (Appendix A). 
Now, for a disordered system, the concept of "bottom of 
the band" has to be taken carefully as the DOS will al- 
ways exhibit Lifshitz exponential tails. In this case, from 
the calculations of the averaged DOS, the bottom of the 
band is found to lie at —At for the paramagnetic (PM) 
(M = 0) case^^ (cfr. Fig. 5). We intend to construct the 
polaronic phase having the paramagnetic, uniform, phase 



and reflects the two competing effects at play: the first is 
the electron's preference for a ferromagnetic background 
accompanied by an energy cost for the localization; the 
second is the reduction of entropy caused by the appear- 
ance of the (fully polarized) magnetic polarons. The last 
contribution, 5cfg, expresses a configurational entropy, 
related to the spacial distribution of the polarons inside 
the system. It is a combinatorial term that can be ap- 
proximated by 




but, since we are working with ^ 1 and the polaronic 
system is below the percolation threshold, it happens to 
be the smallest contribution to AS, allowing us to neglect 
it without important quantitative consequences. 

In writing eq. (6) some important assumptions were 
made regarding the wave function of the electron. Rea- 
soning in terms of the original band state of the electron, 
it is clear that when the magnetic background polarizes, 
the electron energy is lowered by 2t. Thus, the poten- 
tial well is, at most, 2t deep and any bound state will 
always exhibit exponential leaking of the wave function 
to the outside, whereas the energy of the bound state 
in (6) was chosen as the energy of an electron inside an 
infinite potential well. On the other hand, since one can 
think of an effective magnetic coupling as proportional to 
|?/'e(r)|^, the magnetization profile of the polaron should 
also display a smooth variation, whereas in eq. (6) M = 
outside and AI ~ 1 inside the polaron (cfr. Fig. 2). These 
statements amount to say that the exact treatment of 



as reference. Within a virtual-crystal approximation, the 
electron at the bottom of the band has an energy of — 4t, 
and will be extended throughout the system. If a region 
of ferromagnetism develops locally, its reference energy 
in this region will be lowered to — 6i, but there is an ex- 
tra energy that has to be paid if the electron is to become 
confined to this region. For simplicity let us assume that 
the polaron so formed consists of a region inside a cube 
of side R (in units of the lattice parameter, a), inside 
which AI = 1. Obviously, given that in the DEM the 
magnetic interaction is mediated by electron itinerancy, 
one expects that, once the electron localizes inside this 
cube, there will be no magnetism outside, at any temper- 
ature. The variation of free energy per lattice site when 
going from the PM homogeneous phase to this polaronic 
one will be written as 



(6) 

I 

the problem requires a self-consistent calculation of the 
bound state starting from the DEM or perhaps from the 
full Hamiltonian^^'^*. In this sense, eq. (6) is to be un- 
derstood in the spirit of a variational approach, R being 
the variational parameter. 

There are several reasons to expect it to be a good ap- 
proximation: (i) the electron density is very small, mean- 
ing that the overlap of the wave functions of self-trapped 
electrons (and thus the polarons) should be negligible; (ii) 
even if one could devise a full, self-consistent, solution to 
the problem, the exact energy of the bound state is ex- 
pected to differ from the one used in this approximation 
only by numerical factors of 0(1). Specific confirmation 
for this can be found, for instance, in the 1-D calcula- 
tions of Pathak and Satpathy^*, who find that the two 
energies (calculated using the same kind of variational 
wave function used in the current work and the exact 
energy obtained numerically) differ by less than 3 % in 
the limit relevant to our discussion. ; (iii) numerical cal- 
culations strongly favor this approach. To dwell a while 
on this last point, wc studied the effect of a region of 
full polarization embedded in a PM background on the 
electronic spectrum. Performing exact diagonalizations 
of the DEM using spin configurations like the one de- 
picted in Figs. 2 and 3(c), one finds, as a result, the 
appearance of well defined bound states below the con- 
tinuum band. This is clearly seen in Figs. 3(a) and 3(b), 
where the bound states appear with energies and degen- 
eracy coinciding with the value expected for a finite box 

in d-dimensions: ~2c?i cos ^-^^^n^^ At the same 

time, inspecting particular realizations of disorder as the 



AJ^Poi(i?, T) = Atue - 6tne cos 



R+1 



TneR"^ log(25 + 1) - T5cfg(ne, R) 
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(c) 



FIG. 3: Comparison of the disorder- averaged electronic DOS 
for the homogeneous PM phase (hne) with the exact spectrum 
obtained for the same configuration with a polaronic region or 
size R (histogram). Arrows highlight the lowest bound state. 
The results in (a) are for 2D and R — 9, and in (b) for 3D, 
R = 5. The histogram was obtained from the full numerical 
diagonalization of the polaronic configuration. The continu- 
ous line represents the bulk DOS obtained with the recursion 
method. For this comparison we neglected the Berry phase 



in the hopping (i.e. we used at 



,1). Panel (c) shows one 



of those 2D spin configurations in the polaronic phase (bot- 
tom) together with the corresponding exact wavefunction of 
the lowest bound state (top). 



one in Fig. 3(c), one concludes that the wave function 
is clearly localized within the polarized region, thus sup- 
porting our trial function selection for the evaluation of 
the electronic energy. 

The stability of the magnetic polaron is determined by 
the condition AT{Rcq,T) < 0, Rcq being the value that 
minimizes the free energy (6). Important insights can 
be extracted from the direct analytic results obtainable 
in its continuum version (i.e. R 3> 1). In this case the 
equilibrium radius is 



i?eq(T) 



_Tlog(2S'+ 1) 



1/5 



(8) 



increasing at low temperatures as T This power law 
behavior of i?oq with temperature is reminiscent of the 



one found in the context of magnetic semiconductors^^. 
Using this result in the stability condition, one finds the 
stability temperature, Tm, below which the polaronic 
phase appears: 



T„ 



257r3 1og(2S'-f 1) V 5 



t. 



(9) 



These two results convey the essential information for the 
physical situation as one reduces the temperature of our 
system. The high-temperature phase is characterized by 
PM order until T decreases below T,„, at which point, the 
"entropic pressure" is not enough to counteract the en- 
ergy gain and the polaronic phase is stabilized. The tran- 
sition is sharp as a consequence of the one-particle nature 
of the free energy (6), and the polarons set in with a fi- 
nite radius i?eq {Tm)- Continuing the decrease in temper- 
ature, the polaron radius increases until the overlapping 
probability can no longer be ignored, therefrom arising an 
instability towards an homogeneous ferromagnetic (FM) 
phase. An estimate of the Curie temperature, Tq, can 
thus be extracted from a percolation criterion 



yielding 



J r: 



2t:H 
log(25+l) 



5/3 



(10) 



(11) 



where pc is the percolation threshold (for the densities 
we are interested in 3> ri^). For Tc < T < Tm, an 
anomaly in the paramagnetic susceptibility is expected to 
signal the presence of the polarons through an enhanced 
effective moment. 

A natural question can surface at this point regarding 
the fact that, since de Gennes^, we know that, using the 
virtual crystal and the one-particle approach used above, 
a transition between uniform PM and FM phases should 
occur. It is therefore of natural interest to investigate 
how this magnetic transition is altered by the stabiliza- 
tion of the polaronic phase. In order to do that we follow 
the same approach as the one carried by de Gcnnes in 
evaluating the electronic energy as a function of the mag- 
netization, obtaining a description in terms of the free 
energy AJ^pm{T,M). This represents the free-energy 
difference between the homogeneous PM phase and an 
homogeneous FM phase^. 

Since the PM phase is common to A.!Ffm{T,M) and 
AJ-poi {T, M) (6) we can calculate the regions of relative 
stability of the three phases, and draw the phase dia- 
gram shown in Fig. 4. In this plot, T^^^ represents the 
line that would be obtained ignoring the possibility of po- 
laron formation, as de Gennes did; is just the curve 
from eq. (11) corresponding to the percolation criterion. 
Notice that we also included the lines Tc and Tm that 
are the actual transition lines calculated by minimizing 
(6) with respect to the polaron radius and ATpm{T, M) 
with respect to M . Even in such simple phase diagram 



6 



0.010 



o.ooe 



I" 0.004 



0.002 



i63.8 



PM 



Pol 




O.OOOi'- ***^^- — ' 

0.001 



0.002 0.003 0.004 
Electron Density (n ) 



0.005 



FIG. 4: Phase boundaries obtained within the IPM. The 
polaronic phase (Pol) appears as a precursor to the ferromag- 
netism (FM) at low densities (dashed region). The dotted- line 
curve represents Tc obtained in mean field. See the text for 
other notations. The left vertical scale is in units of t, and 
the one on the right in Kelvin, using t — 0.55 eV. 



we observe important physical consequences of the pres- 
ence of the polaronic phase, namely: (a) the polaronic 
phase mediates the transition from an homogeneous PM 
to an homogeneous FM phase, as expected; (b) the tran- 
sition temperature calculated in the mean-field approach 
Tq'^ is notably reduced at low densities by the onset of 
the polarons; (c) the Curie temperature obtained with 
the percolation criterion, Tq , gives a very good estimate 
of Tc (as follows from the superposition of the two curves 
in almost the entire region), supporting the interpreta- 
tion of the ferromagnetic transition arising from polaron 
percolation. 



IV. POLARONIC EVIDENCE IN EUBe 

As mentioned earlier, EuBg is a magnetic metal with 
extremely reduced electron density, and exhibits all the 
characteristic signatures of a polaronic phase in Raman 
scattering measurements^' Such experiments reveal 
that the FM transition at Tc ~ 15 K is preceded by an 
interval of temperatures where the response of the system 
is dominated by the presence of magnetic polarons. 

It has been proposed in Ref. 23 that EuBg is very likely 
a DEM material in the low density regime, characterized 
by a hopping integral of t = 0.55 eV, and a carrier density 
per unit cell rig 3 x 10^"^. Subsequent magneto-optical 
experiments tally with this DE-based interpretation^^. 

If we incorporate such hopping (found to reproduce the 
variation in ujp for EuBg) in the phase diagram for the in- 
dependent polaron model, the result is the one shown in 
Fig. 4, where the absolute temperatures should be read 
in the right vertical axis. It is interesting to compare 
this phase diagram for Ue — 0.003 with the experimental 
evidence from the spin-flip Raman scattering results of 
Ref. 10. The experiment reveals a polaronic region being 
stabilized essentially in the same temperature range as 
the one in Fig. 4 for the appropriate densities. This seems 
to show that our simple description of the polaronic phase 
captures the essential details of the polaron physics in 
EuBg. In particular, we emphasize that, the calculation 
of the Curie temperature based on the de Gennes ap- 
proach (the dotted portion of the straight line in Fig. 4), 
clearly overestimates the actual Tc by a factor of 3 or 
more inside the polaron stability range. As a matter of 
fact, placing ourselves at a density = 0.003 in the di- 
agram of Fig. 4, we find Tc — 17 K. This lies noticeably 
close to the experimental Tq, without adjusting parame- 
ters. 



As discussed above, in the present framework the 
PM — > Pol transition in this diagram is of first order 
and, strictly speaking, the same occurs at the Pol FM 
transition, because the phase boundaries are calculated 
from the relative stability of the two phases. At first sight 
this would mean that, during the Pol — > FM transition, 
the magnetization has a discontinuity at Tc- However, 
from the interpretation of the FM transition in terms 
of polaron percolation, one expects a continuous transi- 
tion, in the sense that the magnetization of the system 
should be weighted by the mass of the infinite percolating 
cluster, which evolves continuously from the percolation 
threshold. Finally, it is interesting to notice the order of 
magnitude of the relevant temperatures and densities for 
this phase. If, for definiteness, one assumes that t ~ 1 cV, 
eqs. (9) and (11) reveal that the stability condition is real- 
ized for densities of ^ 10"'^, and temperatures typically 
in the dozens of Kelvin. As is shown in Appendix A, 
a more realistic approximation for the electronic energy 
doesn't appear to modify these ranges significantly, which 
somehow restrict the range of materials where the effect 
might be realizable. 



V. THE PROBLEM OF PHASE SEPARATION 

The independent polaron model discussed before, is 
based on various approximations that stem from the low 
electronic density of the systems we aim to describe. In 
particular, the same assumption for the electronic en- 
ergies used by de Gennes was employed. One of the 
implications of approximating eq. (4) by Eei{M,ne) ~ 
Ei,{M)ne is that the Curie temperature calculated in 
mean-field for homogeneous PM and FM phases is sim- 
ply proportional to the electronic density. This happens 
because the electronic density, ng, appears together with 
the only energy scale in the problem, t. 

But a more serious question regarding the double ex- 
change in the low density limit is the problem of phase 
separation (PS). It is know that the DEM is unstable 
towards phase separation at low carrier densities, even 
without additional AFM (superexchange) couplings be- 
tween the lattice gpins^^'^^'^^'^'''^^'^^. In order to study 
this aspect of the model in more detail, we will abandon 
the previous single particle approximation for the clcc- 
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FIG. 5: Numerical DOS for the DEM at diflFerent magneti- 
zations, as defined in eq. (21). Tfie arrows signal the position 
of Ef{M — 0) for an electron concentration of Tie — 0.01. 
The inset is a plot of the band edge as a function of the local 
magnetization. 



tronic energy, and calculate this quantity exactly within a 
hybrid thermodynamic approach (HTA) discussed next. 



A. Hybrid Thermodynamic Approach 

The DE Hamiltonian (2) represents a system of classi- 
cal spins that do not interact directly via a conventional 
exchange term. Their interaction comes indirectly from 
the fact that certain configurations of the local spins will 
minimize the energy of the electron gas. At absolute zero 
temperature, it is more or less evident that the ground 
state corresponds to ferromagnetism for any density of 
electrons^*. For other than this specific case, thermo- 
dynamics comes into play. All relevant quantities follow 
from the partition function 



E = P5, Tr 



expl 



'/3{n{S, ) - fiN)) 



(12) 



The integral spans all local spins. Si , and the trace is 
over the electronic degrees of freedom, for which N is 
the number operator. The electrons can be easily traced 
out in this grand canonical ensemble yielding an (exact) 
effective spin Hamiltonian that reads 



(13) 



with 



/3ife// = -^log[l + exp(-/3(i?„(S^O~M))] , (14) 



En{Si) being the eigenenergies of the one-electron 
Hamiltonian. Here, however, lies the origin of the dif- 
ficulties that this system poses to analytical and numer- 
ical approaches. The latter are rather notorious, for one 



might think that once the effective spin Hamiltonian is 
written down, the energy of spin configurations can be 
calculated and the problem can be tackled with usual 
Monte Carlo (MC) techniques. That is indeed so for- 
mally. Unfortunately, the effective Hamiltonian — that 
is, the energy associated with a given spin configuration 
— requires the knowledge of the electronic eigenstates 
associated with such configuration. In other terms, the 
MC methodology would imply the full re-diagonalization 
of the electronic problem at every tentative update of the 
local spin configuration. This is clearly prohibitive, im- 
posing an upper limit of ~ 6^ to 8'^ on the sizes of the 
systems that can be thus studied^^'^^. 

Our approach to this problem tries to circumvent these 
issues through a compromise in which the electronic 
problem is solved exactly and the spin subsystem treated 
within mean-ficld'^^ . It hinges upon the fact that, writing 
Heff in terms of the total electronic DOS 

/oc 
p{E,S.,)log l + cxp{^p{E~fi)) dE, 
-OO 

(15) 

the dependence on the spin configuration is completely 
transferred to the DOS. Since the treatment of the ex- 
act effective Hamiltonian is out of reach, we resort to 
the Bogoliubov-Gibbs inequality for the canonical free 
energy: 



T^-TlogiE)<{H,ff)^^TSt. 



(16) 



Here (...)( means the averages are calculated with a trial 
statistical operator — other than the canonical Maxwell 
distribution — and St is the associated entropy. Since 
our system is expected to exhibit either fcrro or paramag- 
netism, the simplest suitable choice is the one generated 
by the uniform mean-field Hamiltonian 



N 



(17) 



where /i is a variational parameter used to minimize the 
inequality (16). Since all averages are now done with 
regard to this TYj, we have: 



M = {^S, j^ ^ J VS, exp{-pnt)S^ 

- (coth(/3/.) - ^) 
= C{Ph) Uz , 
where C{x) is the familiar Langcvin function, 

sinh(/i) 



TSt = log 



hM , 



(18) 



(19) 



and 



p{E,S^)j log[l+exp(-/3(S-M))Jrf£;. 

(20) 
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Therefore, the trouble in calculating the equilibrium free 
energy now boils down to the computation of {p{E, Si ))^. 
We know^^'™'^^ that the recursive method can be used to 
obtain the exact p{E, Si ) for a given configuration {Si }. 
It is then a matter of straightforward statistics to ob- 
tain the averaged DOS. Hence, the electronic problem is 
still treated exactly for every configuration of local spins. 
These configurations are generated with the probability 
distribution ~ cxp(— /3ft, S'i ^) and, since /3/i and A4 
arc univocally related through (18), we simplify the no- 
tation and write 

{p{Ej,))^=p{E,M) (21) 

whenever we refer to the DOS averaged over configura- 
tions of disorder compatible with an average magnetiza- 
tion M. 

We also recall that when Ep T, the logarithm in 
Eq. (20) can be replaced by 0(^ — E) meaning that the 
thermal excitations of the electronic subsystem can be 
neglected to a great extent. According to our earlier dis- 
cussion, we are interested in cases for which t 1 cV and 
Tq ^ t. This suggests a "zero-temperature" description 
of the electronic system, the thcrmal/entropic effects be- 
ing assigned entirely to the spin subsystem^^. 



B. Canonical Free Energy and PhaseDiagram 

Unless whenever stated otherwise, throughout this sec- 
tion we will be concerned only with the homogeneous 
phases (PM and FM) of the DEM. In the HTA one tries 
to trace the electrons out of the problem and obtain an 
effective Hamiltonian for the lattice spins, in such a way 
that the partition function becomes simply (14): 

S = j VS,e^p{-l3H,ff{S,)) . 

This effective Hamiltonian contributes to the total free 
energy through eqs. (16), (19) and (20), and we have, at 
the end: 

T[M, ne) = / ep(e, M)dt - TS(M) . (22) 

The free energy in eq. (22) should be minimized with re- 
spect to M to obtain the equilibrium free energy !Feq{ne)- 
We work at constant electron density and our focus in 
the cases where T -^t is implicitly assumed in the omis- 
sion of the electronic entropy. The phase diagram that 
emerges from the minimization of (22) is drawn in Fig. 6 
and reproduccs^^ the results obtained by Ref. 35. The 
PM-FM transition happens to be continuous for all den- 
sities and the deviations from the result obtained with 
the bottom of the band, de Gennes-like, approximation 
are quite evident. It turns out, however, that this phase 
diagram is incomplete. 




"e 

FIG. 6: Phase diagram of the DEM obtained within 
the HTA. For comparison, the result for Tc obtained by 
de Gennes is also plotted (dashed line) [cfr. Ref. 35]. 



C. The Essence of the Problem 

Reflecting the relative stability of the two homoge- 
neous phases, the plot in Fig. 6 says nothing regarding 
density fluctuations. To determine whether the system 
is unstable in relation to density fluctuations, we need 
to look at the behavior of the equilibrium free energy, 
-^eij('^e), namely its dependence on the electronic den- 
sity. Thermodynamic stability requires J-eq(ne) to be a 
globally upward convex function, so that the compress- 
ibility, K, = n~^dn/dp,, is never negative. Whenever this 
condition is violated and the density of the total system 
is kept constant, it will naturally segregate into two dis- 
tinct phases in such a way that guarantees the restoration 
of convexity in the resulting free energy. The equilibrium 
state in this phase-separated region can be obtained by 
the so-called Maxwell construction, which, geometrically, 
is tantamount to substituting the underlying J^eq{ne) by 
the envelope of all inferior tangents^^. A sketch of the 
situation is presented in Fig. 7. 

It so happens that, when the behavior of the isother- 
mals of J^eq{ne) for the DEM calculated using eq. (22) 




FIG. 7: Free energy for a system unstable towards den- 
sity fluctuations. The system is unstable for ua < n < ub 
and segregates between regions of density ha and ub- The 
Maxwell construction restoring the curvature of J-eqijie) is 
represented by the dashed lines. 
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dashed (red) curve is the chemical potential resulting 
from the free energy after the Maxwell construction. Be- 
low n__ and above n+ the result is the same, but in be- 
tween it is simply an horizontal line. Altough not shown 
in the figure, the position of this horizontal part is such 
that the areas Ai and A2 are exactly equal, as it should 
happen^^. The third curve (circles) shows the chemical 
potential obtained my minimizing the free energy in the 
grand canonical ensemble. In this case, since not n but 
/i is kept constant, the phase separation instability arises 
naturally through a discontinuity in the ^{n) curve. The 
discontinuity appears precisely in the region of densities 
where the Maxwell construction is in effect, and is just 
another way to observe this phenomenon. 



0.01 0.02 0.03 0.04 0.05 0.01 0.02 0.03 0.04 0.05 

Density - n > 



FIG. 8: Details of the Maxwell construction for the DEM 
at T = O.OSt . The result of the construction is barely dis- 
cernible from the underlying J-eqijie) (I) and, for clarity their 
difference is plotted on panel II. For comparison, the chemi- 
cal potential extracted from the results in I is shown together 
with calculations in the grand canonical ensemble (III). Panel 
IV shows d^Tfiq/dni . 



is analyzed, the signature of this instability emerges at 
low densities by the violation of the global convexity. A 
typical result is shown in detail in Fig. 8. Since there 
is a considerable amount of information in this figure let 
us go through it in detail (everything is calculated at a 
constant temperature, T = 0.03i). 

In the first panel (I) we show the equilibrium free en- 
ergy J-eq{ne) as a function of the electronic density 
(continuous/black line). Despite looking like a straight 
line, it does have a slight upwards curvature at the low- 
est densities, downwards at intermediate densities, and 
upwards again at higher densities. This can be seen in 
the curve d^Teq/dn^, plotted in panel (IV). Thus, we 
have an instability and a Maxwell construction has to 
be done in order to find the true equilibrium free en- 
ergy of the system. The Maxwell construction is the 
dashed/red line in panel (I). Since the effect is rather 
subtle, we plot the difference between J^eqine) and the 
free energy after the Maxwell construction in panel (II). 
In panel (I) we also show the density corresponding to 
the PM-FM transition at this temperature, signaled by 
the dot-dashed vertical line (cfr. Fig. 6). If we call this 
particular density ne{Tc = 0.03), then it is clear that 
n_ < rif-iTc = 0.03) < n_|_. Then, although the rela- 
tive difference between the free energies is only ~ 1%, 
it is qualitatively significant because the system exibiths 
coexistence of PM and FM, each of the coexisting ther- 
modynamic phases having its own electron density. 

In panel (III) we present the chemical potential 
calculated in three different ways. The first one 
(black/continuous) is simply the curve corresponding to 
dFeq{n)/dn, and the instability is also clearly seen here 
since /x(n) should be monotonously increasing. The 
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FIG. 9: (a) Isothermals in the ^ — rie plane for selected tem- 
peratures. The coexistence curves n_ (squares) and n+ (cir- 
cles) were obtained from the Maxwell construction, whereas 
the fi{ne) curves are calculated in the grand canonical ensem- 
ble, for comparison, (b) Detail of the low density region of the 
updated T — phase diagram. The phase separated region 
(PS) is bounded by the curves n-{T) and n+{T) obtained 
from the Maxwell construction. 
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Actually it is quite instructive to pursue in some more 
detail this complementarity between the canonical and 
grand canonical treatments for this particular case. Our 
scenario is completely analogous to the well known be- 
havior of the van der Waals isothermals for the liquid- 
gas transition^^. This is best understood with reference 
to the plots of /i(ne) in Fig- 9(a) at several temperatures. 
Just as in the P — v phase diagram of the gas, the fi{ne) 
curve is monotonous for temperatures above some crit- 
ical value. Below this point, the compressibility of the 
system {drie/d^) becomes negative in some density in- 
terval implying an instability. The coexistence curve is 
then defined by n_ and obtained from the Maxwell 
construction and, naturally, coincides with the jump in 
the density obtained in the grand canonical calculation. 
Despite the analogy and resemblance of the coexistence 
curve for the DEM and the van der Waals gas, there is 
an important qualitative detail not present in the latter 
case: the coexistence curve for the DEM is re-entrant. 

The stability analysis was performed for all the tem- 
peratures and densities shown in the phase diagram of 
Fig. 6. The system is unstable towards phase separation 
at low densities, and the qualitative behavior of the ther- 
modynamic functions is always the one discussed above. 
As a consequence, one obtains the updated phase dia- 
gram presented in Fig. 9(b). 

The information conveyed by this diagram can be 
translated as follows: for densities above ~ 0.04, the 
system is homogeneous and exhibits a magnetic transi- 
tion at Tc{n^). If the density is lower than 0.01 
the system is homogeneous and PM at high tempera- 
tures until n-{T) is reached. At that point the homo- 
geneous phase is no longer sustainable and two phases 
start to coexist: one with density ?^-(T) and PM to- 
gether with another of density n+(T) and FM. Below 
a temperature signaled as Tps-i, n^{T) = and thus 
the PM portion of the system is devoid of electrons^^. If 
0.01 ^ rie < 0.04, the re-entrant nature of the coexis- 
tence curve means that the system can become an homo- 
geneous ferromagnet below Tc{ne) but still segregate at 
lower temperatures whenever the condition rig = n^{T) 
is met. In any case, the magnetization of the system is 
always continuous - a simple consequence of the relation 
between the densities and volume fractions of each phase: 



at some high temperature T > Ti and has a given den- 
sity n. Under these circumstances, the phase diagram 
states that the equilibrium corresponds to an homoge- 
neous PM phase. We can lower the temperature until 
Ti is reached, at which point an instability arises. Ex- 
actly at r = Ti there will be a segregation between a 
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PS 







«, n-n. 



FIG. 10: Schematic representation of Fig. 9(b). 



PM phase with density ni_ = n and another, FM, with 
density In order to satisfy the constraint (23), the 

volume fraction of the FM phase will be Vi+ = 0, at 
T = Ti . A slight decrease in the temperature from Ti to 
T2 will reorganize the system so that the PM phase with 
density n2- now coexists with a FM phase of density 712+ . 
For T2 very close to Ti, the density of the PM phase is 
just slightly different from the global density: 



n2- 



n — e , 



(24) 



where e is a small quantity. This determines the volume 
fractions to be 



V2+ 



"-2+ 



(25) 



Evidently, this means that most of the electrons still 
remain in the PM phase, and just a few populate the FM 
regions. But this poses a problem. The Coulomb interac- 
tion will certainly prevent the accumulation of charges in 
a very small volume and the system must remain neutral. 



V- 



V 



(23) 



This constraint introduces some peculiarities regarding 
the nature of the phase separated state. To understand 
that, in Fig. 10 we draw a sketch of the phase separated 
region in the phase diagram of Fig. 9(b). Using this 
sketch as reference, assume that our system is initially 



D. Electrostatic Suppression of Phase Separation 

The Maxwell construction is inexpressive with regards 
to the spatial organization of the phase separated state. 
This follows from the fact that the Maxwell construction 
for two coexisting phases of densities n+ and n_ amounts 
formally to saying that. 



J 



V+ V- 



- x) , 



(26) 
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FIG. 11: Depiction of the Wigner-Seitz construction dis- 
cussed in the text. Each cell is replaced by an effective spher- 
ical cell characterized by the same densities. 



where x ~ V+/V and V-/V represents the volume frac- 
tion of the two coexisting phases. This is simply a 
hncar interpolation between the two values F(n^) and 
F(n_), as illustrated in Fig. 7. Therefore, the result- 
ing free energy, J-Maxweii, corresponds to the situation 
where we have a system made of two independent ther- 
modynamic components. In particular, the Maxwell pre- 
scription above says nothing about the way the system 
reorganizes when it phase-separates. This can only come 
from additional interaction terms that should be added 
to the right hand side of eq. (26), in order to include in- 
teraction, surface, boundary and perhaps other relevant 
effects or correlations. In our case the two phases have 
different electronic densities, both different from the ho- 
mogeneous density, the latter satisfying 



n+x + n_ (1 — x) 



(27) 



Obviously, with mobile negative charges, the system can 
indeed adjust itself by accumulating electrons in some 
regions, and depleting them from others. But since the 
background of positive atomic charges stays essentially 
immutable and homogeneous, this means that the total 
charge of each phase is not zero. Coulomb interactions 
are therefore crucial. Under this assumption we now in- 
vestigate two important corrections for the free energy in 
the phase separated regime. 



1. Electrostatic Correction 

Given that we are not assuming any sort of anisotropy, 
the electrostatic constraint will most certainly favor the 
development of bubbles of the FM, high density phase, 
dispersed in the PM, low density phase. Such scenario 
is schematically depicted in Fig. 11. The charge den- 
sity is assumed uniform and continuous inside each FM 
bubble and across the PM background. To calculate the 
electrostatic energy associated with this charge distribu- 
tion wc take notice to the fact that, on the grounds of 
overall charge neutrality, it should be possible to find an 
appropriate neighborhood around each FM bubble such 
that the total charge on the bubble plus PM neighbor- 
hood adds to zero. In Fig. 11 this is represented by the 



wavy lines that, in this way, define cells of charge neu- 
trality. Following Wigner^^, these cells are replaced by 
the equivalent Wigner-Seitz (WS) spherical cell contain- 
ing the same volume fractions of FM and PM phases 
(as shown on the right-hand side of the diagram), which 
means that 



Ri 



xRt 



(28) 



For each WS cell, the total electrostatic energy is cal- 
culated considering three terms 



Ec = U. 



U+-. 



(29) 



with the first accounting for the electrostatic self-energy 
of the region, the second the self-energy of the "-" 
region, and the last the mutual electrostatic interaction 
between the two. All of them are calculated within clas- 
sical electrostatic theory assuming uniform charge distri- 
butions in the two regions. Hence 



3Ql 

5R+ 



(30) 



and represents the electrostatic energy of the inner sphere 
containing the "-f " region; 



3Q2 



Rt - Ri 



Rt - i?3 



Rl 



Rl-Rl 



(31) 



stands for the energy of the outer shell of the "-" region, 
and 



(32) 



is the electrostatic interaction between the inner sphere 
and outer shell. In the above Q± are the total charges 
(positive background -I- electrons) inside the two regions, 
and R± the respective radii, as depicted in Fig. 11. For 
reasons regarding numerical stability, in the following we 
will take always n_ = 0. Looking at Fig. 9(b), this comes 
as a natural approximation because n+ ^ n_ for Tpsi < 
T < Tps2, and becomes exact for T < Tpsi- Hence, 
using the identities 

Q+ — e{n — 77.+ )V+ = en{x — 1)V 

V. = -jRi 

Q_ = e{n - 7i_)F_ = en{l - x)V 
V_^^{R^_-Rl) 



(33) 



the Coulomb term can be cast as 

Ea^-eWRl , 

which yields an energy per unit of volume of 

Ec 2 2 2 p2 2-_3xV3 + x 

ec = -rr- = - e n nR 

V 5 X 



(34) 



(35) 
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For simplicity we replaced by R and keep this lighter 
notation below. The bubble radius, R, is a variational 
parameter. The result (35) implies that, for a given vol- 
ume fraction of the two phases, say x, the electrostati- 
cally favorable situation is to shrink the FM bubbles to 
an arbitrarily small size, meaning that R 0^^. But this 
treatment is yet incomplete, inasmuch as the considera- 
tion of finite-sized bubbles of electrons requires another 
correction, of different nature, to the ground-state en- 
ergy. 



correction to the electronic DOS of the free electron gas 
which, as discussed in Appendix B, leads to a correction 
to the ground state energy per electron reading 



E{R,n+) E{oo,n+) 



15 

^ 16 V 671 



1/3 



1 

'r 



(36) 



The term in R~^ is the correction for each bubble. To ob- 
tain the total correction we just multiply by the number 
of WS cells, obtaining the total energy per unit volume 



2. Phase Space Correction 

The electronic contribution to the free energy (22) is 
calculated in the thermodynamic limit. In the phase sep- 
arated regime, the electron rich bubbles are expected to 
be of relatively small size. Therefore, one cannot rely on 
the electronic energy calculated in the thermodynamic 
limit and the need to introduce finite size corrections 
arises. 

The leading correction to the energy of an electron 
gas confined to a finite sized volume comes through the 



E{R,n^,x) E{oo,n+) 



V 



V 



-{-) 
16 \QnJ 



1/3 2;l/3 



R 



(37) 

As expected, the phase space correction acts to the effect 
of rising the ground state energy of the electron gas. This 
induces a tendency opposed to the one embodied in the 
electrostatic term (35), and the two should balance at 
some optimum value of R. 

The free energy in the phase separated regime (26) 
needs to be updated for these two corrections that go 
beyond the Maxwell construction. The result is 



J 



J^Maxweii{n) = T{n+)x + F{n-){\ - a;) -|- 



V 16 V6ny R 



/3 2;4/3 



(38) 



Since the FM radius is entering only in the correcting 
terms, the minimization with respect to R can be per- 



formed at once, and the final result is 



J 



J^Maxweiiin) = T{n+)x + J'{n^){l - x) + 



3 

22/3 



2e2 
5 a 



1/3 



5 ^ ' 16 



/ \ 1/3' 
— ) 


2/3 


'2 + a; -3x^/3' 


1/3 


\6neJ 




^5/3 





(39) 



This last expression is the free energy per site of the 
original lattice when phase separation is in effect, a is 
the lattice parameter, rie the electron density per unit cell 
of the crystal, and t is the hopping integral. Naturally, 
when there is no PS instability, x = 1 by definition and 
the above reduces to J-{ne) as one certainly expects. In 
the PS regime, the equilibrium radius of the electron rich 
FM bubbles satisfies 



4 „^ 15 at 
- 7rri+i?3 = 

3 + 16 e2 



6n rif 



1/3 



1 



2 + x- 3a;i/3 ' 



(40) 



and this relation can be used, for instance, to inspect the 
typical number of electrons inside each FM bubble. 



E. Consequences for the Phase Diagram 



The natural question is now: what happens when the 
equilibrium free energy is recalculated with these correc- 
tions? Namely, we want to know wheather the PS insta- 
bility persists when the Maxwell construction is updated 
according to (39). It is useful to have a tuning param- 
eter that interpolates between the case in (39) and the 
previous calculation where electrostatic and localization 
effects were disregarded. With that purpose, we intro- 
duce the dielectric constant, e^, that renormalizes the 
electron charge as e^/sr in the expressions above. 

By varying Er between 1 and 1000 the curves in Fig. 12 
were obtained. In this plot, we are focusing on the low 





FIG. 12: Phase diagram calculated within the HTA with 
a corrected Maxwell construction according to eq. (38). The 
arrow marks the polaronic stability region, represented by the 
dashed area at low density and temperature. Notice how the 
region of phase separation lies inside the polaronic bubble for 
dielectric constants of £r ~ 1. 



density region of the phase diagram where PS occurs. 
The red (circles) curve pertains to the case £^—'■00 (or, 
equivalently, zero electronic charge) and is again the re- 
sult shown before in the phase diagram of Fig. 9(b)^^. 
The figure is transparent as to what happens when the 
Coulomb interaction is turned on: the PS region is pro- 
gressively reduced! Not only that but it is clear that 
an overwhelming shrinking of the PS region takes place 
when e = 1, which is a reasonable value, considering the 
low electronic densities. 

Thus, the consideration of the free energy (39), cor- 
rected for the effects arising as a consequence of charge 
accumulation, leads to the suppression of the PS instabil- 
ity. The electrostatic payoff involved in the segregation 
leads the system to phase-separate only at much lower 
densities and/or temperatures. Just how low these are is 
controlled by the effective electronic charge. 



F. Phase Separation and Magnetic Polarons 

There is a question of relevance that we have been post- 
poning since the beginning of this section on the problem 
of PS. In Sec. IV calculated the phase diagram of Fig. 4 
describing the stability conditions for free magnetic po- 
larons in the DEM. From the polaronic phase diagram 
follows that magnetic polarons are only stable at con- 
siderably low temperatures and densities. In fact much 
lower than the temperatures and densities at which the 
PS instability sets in. The reader might have noticed 
since Fig. 10 that the density and temperature scales for 
the PS bubble are much higher than the scales for the 
polaronic bubble. More precisely, the polaronic phase 
lies well inside the PS region when the corrections to 
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FIG. 13: The number of carriers inside each FM bubble 
calculated according to (40) for selected temperatures and 
dielectric constants. 

the Maxwell construction are ignored. These two regions 
can be seen in perspective in Fig. 12, where the polaronic 
stability region is highlighted by the dashed region in the 
lower left corner, and completely inside the PS region for 
e ^ 00. 

This is clearly a problem to our arguments concerning 
the polaronic phase. The polaronic stability has been de- 
termined by studying its relative stability with regards to 
an homogeneous FM phase. The diagram above is say- 
ing that, at such low densities, there are no homogeneous 
phases — the system phase separates! So the study of 
polaronic stability has a problem. 

But this is only if the electrostatic effects are ignored. 
With their inclusion, the PS region retreats to lower and 
lower densities and, as the figure documents, for Sr ~ 1, 
it rests already completely inside the polaronic region. 
So, it seems that the problems above with the polaronic 
phase have just diminished. 

The connection between PS and magnetic polarons is 
indeed remarkably close. The ferromagnetic droplets, as- 
sociated with a localization energy for the electrons in a 
restricted volume, arc nothing more than our descrip- 
tion of the magnetic polaron. So, in this sense, the PS 
regime studied here and the magnetic polarons are dif- 
ferent perspectives of the same physical concept. One 
of the differences is that, while the magnetic polaron is 
defined as a FM droplet with a single electron, the PS 
regime allows for droplets with many electrons. 

To explore this further, it is interesting to know the 
number of electrons inside each FM bubble in the PS 
regime. Some typical results are plotted in Fig. 13 for 
the same used before. The most remarkable fact about 
these curves is that the one pertaining to = 1 is of the 
order of unity. Therefore, there is essentially one electron 
per FM droplet. In addition, the temperatures at which 
PS occurs for this value of Sr are so low that the FM 
droplets are very nearly full polarization. But, a fully 
polarized FM droplet with one electron inside is just our 
definition of magnetic polaron! Thus, the peculiarities 
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of the phase segregation in this system are completely 
consistent with the magnetic polaron picture and, with 
that respect, the similarity between the shape of the two 
phase diagrams (Figs. 9(b) and 16(a)) seems hardly co- 
incidental. 



G. General Argument Regarding Phase Separation 



Although wc have been focusing so far on the specifici- 
ties of the phase separation instability in the DEM, phase 
separation is quite a general phenomenon in thermody- 
namics. To conclude our discussion we put forward an 
argument showing that, if the contribution of the entropy 
of the electronic gas can be neglected in the total free 
energy, then phase separation is ubiquitous in electronic 
systems whose bandwidth is magnetization dependent. To 
see how this comes about, recall that the free energy (22) 
is given by only two contributions: the electronic ground 
state energy plus an entropic term attributed uniquely 
to the local moments. For reasons that will be clear in a 
moment, this argument unfolds more clearly if we work in 
the grand canonical ensemble, where the electronic chem- 
ical potential, is held constant. In this case the free 
energy (22) is simply replaced by the grand canonical po- 
tential, and is akin to replacing Ep — > fi and e e — fj, 
there. 

Consider now the following facts, regarding the rele- 
vant thermodynamic quantities seen as functions of the 
magnetization: 

1. The magnetic entropy (19) is monotonous and 
downward convex, for all the domain of A/; 

2. For a given chemical potential, ^, the electronic en- 
ergy is monotonous and downward convex through- 
out the entire domain of variation of Af ; 

3. The electronic density is monotonous and upward 
convex. 

Statement 1 follows from the fact that S{M) is a proper 
thermodynamic entropy for a magnetic system, having 
all the required analytical properties. Point 3 is a trivial 
consequence of the electronic density being the integrated 
DOS. Point 2 can be understood from the fact that the 
electronic bandwidth is monotonous with M (cfr. Fig. 5). 
There is a subtlety however in that, if fi happens to be be- 
low the band edge at M ~ 0, then, the electronic energy 
will be identically zero until some critical magnetization, 
say M*, is reached for which the band edge coincides 
with /i. For higher magnetizations, the energy decreases 
and the overall shape is as depicted in Fig. 14. Since 
the same plateau is present in the electronic density, for 
exactly the same reasons, the result for the grand canoni- 
cal potential will be something like the solid curve drawn 
schematically in the bottom frame panel of Fig. 14. Ev- 
idently, there will be a temperature at which the min- 
imum of this curve exactly touches the horizontal axis 




M* ^ M 1 M 



Eel \ 

FIG. 14: Schematic variation of DOS, entropy and electronic 
energy with M at constant fi. In the top panel we illustrate 
the changes in the DOS with increasing magnetization, and 
the corresponding increase in bandwidth. In the lower panel 
we illustrate the variation of the electronic energy and mag- 
netic entropy, and how they combine to yield the free energy 
depicted by the black (solid) line. 

(as depicted), thus precipitating a first order transition. 
At the precise temperature, T, at which this happens, 
the system stays undecided as to which state it should 
have because the thermodynamic potentials at Af = 
and M = Mps are degenerate. Since /i is kept constant, 
M = and M = Mps correspond to different electronic 
densities. This is but our phase separation instability 
seen from the grand canonical perspective. 

The important thing here is that nothing in this ar- 
gument mentions the details of the specific model under 
consideration, and therefore is valid as long as the basic 
assumptions remain valid. In particular, the magneti- 
zation is as good as any other suitable thermodynamic 
parameter, and, thus, the arguments extends to any ap- 
propriate classical variable coupled to the electronic en- 
ergy as the magnetization is in our specific case. 

That phase separation has to emerge always follows 
from the fact that one can always place /i in between the 
band edges at Af = and M = 1. So, if temperatures are 
so low that the fermionic entropy can be disregarded, this 
kind of treatment should always yield a phase separated 
regime at low densities. 



VI. THE REGIME OF ULTRA LOW DENSITIES 

Having clarified the issue of phase separation and its 
connections with the polaronic phase, we return now to 
the polaronic description. The consideration of electron- 
electron interaction, as done in the previous section, 
amounts to effectively describing the electrons in terms 
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of the Hamiltonian 



Hde = - ^ O'ijdldj + h. c. 

1>J ' ' 



(41) 



where the presence of the Coulomb term is now exphcit. 
The presence of the full Coulomb interaction in Eq. (41) 
is to be understood in connection with the case of ex- 
tremely reduced electron density. In a conventional Fermi 
liquid the high density of electrons makes the screen- 
ing process very effective, and the effect of the electron- 
electron interactions can be absorbed into the renormal- 
ization of physical quantities such as the effective mass. 
In a very diluted electron gas the Coulomb interaction 
cannot be addressed meaningfully in this way. In the 
language of the one component plasma, this can be un- 
derstood with reference to the dimensionless parameter 
i~s = fo/o-Oy where — {S/AnpY^^ {p being the volumet- 
ric density of electrons) is the average distance between 
electrons, and ao is the Bohr radius. While the kinetic en- 
ergy scales as l/r^, the potential energy varies as l/r^, 
and; therefore dominates in the low density (rs — > oo) 
regime^^. 

It is well known since Wigner^" that, under such cir- 
cumstances, the electrons arrange themselves in a reg- 
ular lattice, and electron localization occurs. If a so- 
called Wigner crystal is realizable in a system described 
by Eq. (41) then, as in any solid, there is zero point mo- 
tion of the electrons about their equilibrium positions. 
If the electron itinerates among several unit cells dur- 
ing this zero point motion, the magnetic coupling, Jh, 
leads to the local polarization of the lattice spins, pro- 
ducing a bound magnetic polaron and generating a po- 
laronic Wigner crystal. This process is obviously limited 
by the melting of the electronic solid and, thus, the ques- 
tion arises of how to describe these two tendencies for 
the polaron formation and the interplay of Coulomb and 
magnetic interactions. We present estimates regarding 
the stability of this polaronic Wigner crystal below. 



A. Polaronic Wigner Crystal 

Following the Wigner-Seitz approximation^^, the 
Wigner crystal unit cell (much larger than the original 
lattice spacing, a) is approximated by an electrically neu- 
tral spherical volume, inside which the ionic charge den- 
sity is homogeneous. The electrostatic potential energy 
then depends only upon r: the distance of the electron 
from the center of the cell. The Hamiltonian for an elec- 
tron in this uniform charge and spin background is then: 

nw^-Qt-'^'- + ^ + \m{u;l + u^)r^ , (42) 

where p is the electron momentum, m — 1/(20^1) the 
effective electron mass, and = w^/S = /mr^, where 
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FIG. 15: Phase diagram of the diluted DEM model as a 
function of temperature and density {t — 0.5 eV, a — 4 A): 
paramagnetic (PM), polaronic (Pol), polaronic Wigner solid 
(WS), and ferromagnetic (FM). Notice that the axes are pre- 
sented on a logarithmic scale (cfr. Fig. 4). 



LUp is the plasma frequency. In Eq. (42) to is the frequency 
of the confining potential due to the DEM mechanism, 
and is a variational parameter. The radius of the mag- 
netic polaron in the Wigner crystal relates to to by: 



R = 



(43) 



where 51 = y^uu'^ + ui'^ is the total frequency of oscillation 
of the electron. Notice that the ground state energy of 
Eq. (42) is Eo = -6t - 3e'^/{2ro) + 3n/2 and hence the 
relative gain in free energy is: 



AJ^VFP = -2t + 3 {n - uje) /2 + 4TTR^T\og{2S + 1) /3 
Minimization of Eq. (44) with respect to uj gives 
Rs {T*/Tf'^ ,T>T* 



R{T) 



Rs 



(44) 
(45) 



where Rs = [3t/ {uje)]^^^ is the saturation radius, and 
T*/t = 9/[47ri?|log(2S'+ 1)] is the temperature below 
which the polaron radius saturates due to the interplay 
between the DEM and the Coulomb interaction. 



B. Wigner Crystal Melting 

It is clear that the previous results are valid for temper- 
atures so low as not to melt the Wigner solid. The calcu- 
lations for the independent polaron model reveal that the 
temperatures for polaron stability are already typically 
small for reasonable values of t, but the electronic solid 
is much more sensitive to the temperature. The Wigner 
crystal melting temperature, Tm, can be estimated from 
the Lindemann's criteria ''"'^: Tm ~ 0.01{e^ /a)nl^^. It is 
known from several numerical calculations''^'''^'^^ on the 
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stability of the one component plasma, that the max- 
imum densities and temperatures at which the Wigner 
crystal can exist correspond to ~ 50-100, and T ^ 
10 K. Values of rs ^ 50-100 correspond to rig ^ 
10~6-10"^ for t = leV and a = 4A. The region of 
stability of the polaronic crystal is show in Fig. 15. Due 
to the absence of magnetic interactions between different 
polarons, the Wigner crystal should be a superparamag- 
net: the local moments within the zero point radius of 
the electrons are expected to respond collectively. In the 
presence of other long-range interactions (such as dipole- 
dipole) the polaronic Wigner crystal can exhibit long 
range magnetic order. Increasing the electron density 
at T = causes the Wigner crystal to quantum melt at a 
critical density with two possible outcomes: a paramag- 
netic polaronic Fermi liquid or a fully polarized ferromag- 
net. In both cases the carriers are mobile and can screen 
the long-range part of the Coulomb interaction leading 
to a Fermi liquid state. At finite temperatures, where the 
electron state cannot be described by the zero point mo- 
tions implicit in Eq. (42) alone, the crystal should follow 
the features of the phase diagram for the electron gas^^ . 
The characterization of the system in the neighborhood 
of the melting point, where the presence of a polaron 
liquid is plausible [Fig. 15], is restrictively hard, even 
for the simple electron gas. Far from this region, where 
the electron density is high enough to make the screen- 
ing process effective, one expects to retrieve the behavior 
obtained before within the independent polaron model, 
and discussed previously. 

VII. CONCLUSIONS 

The pure double exchange model Hamiltonian (2) dis- 
plays a rich phase diagram, even at the lowest densi- 
ties. This is in itself noteworthy insofar as we did not 
introduce any additional competing interactions, such as 
direct magnetic couplings among the local spins. In ad- 
dition to its intrinsic tendency for a FM transition at 
low temperatures, the low density region of the phase 
diagram displays an instability towards phase separa- 
tion. This phase separated regime is characterized by 
FM droplets rich in electrons, embedded in a PM back- 
ground, essentially depleted of electrons. The consider- 
ation of electron-electron interactions shows that, as ex- 
pected, the phase separated regime is highly suppressed 
for meaningful! magnitudes of the electrostatic correc- 
tion. This suppression is such that the effective number 
of electrons inside each FM droplet becomes of the or- 
der of unity. Given that, by construction, our WS cells 
containing the FM droplet are neutral, this situation of 
nearly one electron per droplet can be alternatively ad- 
dressed in terms of non-interacting magnetic polarons. 

This justifies our initial approach to address the pola- 
ronic phase in the low density double exchange model. 
Within the independent polaron model developed in 
Sec. Ill, below a critical density the PM-FM transition is 



mediated by a polaronic phase (cfr. Fig. 4). One conse- 
quence of this is that the Curie temperature of the system 
is much lower than one would obtain if based only on a 
PM-FM mean-field approximation. Moreover, by consid- 
ering the model parameters already used in Ref. 23 to de- 
scribe other properties of EuBg from a DE point of view, 
we obtain the range of temperatures Tm < T < Tq for the 
presence of the polaronic phase in agreement with exper- 
imental observations, and without additional adjustable 
parameters. 

Further down in the density scales we estimated the 
conditions for the stability of a polaronic Wigner crys- 
tal. This phase seems plausible in the ultra diluted sit- 
uation where Wigner cristalization of the electron gas is 
expected. In this case the zero point motion of the elec- 
trons can still provide enough itinerancy to polarize the 
neighboring local moments, generating a crystal of bound 
magnetic polarons. Nevertheless this regime, although 
rather appealing from the theoretical point of view, is 
certainly difficult to reach experimentally on account of 
the reduced temperature and density scales involved. 

The current results complement the ones in Ref. 23 
that pertain, mostly, to the evolution of the system once 
the homogeneous FM phases sets in, whereas now we 
have addressed how the transition from homogeneous PM 
to the onset of homogeneous FM takes place. We can 
then interpret the ferromagnetic transition in EuBg as 
being precipitated by the merging of magnetic polarons 
which attain the percolation threshold close to Tq- At 
the same time, these results lend additional support to 
the interpretation of the phenomenology of EuBg from a 
double exchange perspective. 
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APPENDIX A: EFFECTS OF FINITE BAND 
FILLING ON POLARON STABILITY 

In the main discussion of the polaronic physics in the 
DEM, the simplification is made of considering that the 
electronic energy is simply accounted by the energy at the 
bottom of the band, multiplied by the electron density 

(5). 

In order to account for the finite electronic density, we 
rely on a rigid band approximation for the DOS. This 
means that we calculate the electronic energy for a given 
density of polarons, rip, assuming the DOS in the conduc- 
tion band doesn't change significantly^^. The free energy 
per lattice site then becomes 
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Fpoi{R, M, Up) = —Qtup cos 



R+l 



Ep{nc—np,M) 



BUM) 



p{e, M)ede + (1 - npR^)TS{M, S) . 



r 



In this approximation, the electronic energy is counted 
essentially by transferring electrons from the conduction 
band to the bound states. The Fermi energy satisfies 
Epirie — Up.M) > Ei,{M), reflecting the existence of 
TT-e — Hp electrons in the band. An important differ- 
ence relative to the case of the empty band considered 
in Sec. Ill follows: since we are allowing the existence of 
ne — Up states extended throughout the system, the non- 
polaronic part can still be ferromagnetic below some tem- 
perature. So, in principle, the magnetic polarons could 
be embedded in a background with finite spin polariza- 
tion, and, therefore, we have to introduce the magneti- 
zation of the background, M, as a third variational pa- 
rameter, together with R and Up. The minimization of 
Al with respect to these parameters, and the comparison 
of the resulting equilibrium free energy with the homo- 
geneous case (22), produces the phase diagram displayed 
in Fig. 16(a). The overall qualitative features obtained 
previously in Fig. 4 remain basically the same, the im- 
portant differences now being: (i) The PM-FM transi- 
tion curve (dashed line) now places Tc for homogeneous 
phases at higher temperatures than the ones obtained 
with the de Gennes treatment; (ii) The polaron stabil- 
ity temperature, T^, is seen to increase with density if 
Ue ^ 0.001, just as expected because, the higher the den- 
sity, the higher the Fermi energy and the more favor- 
able it becomes to create a polaron for the same price 
in entropy; (iii) The reentrancc of the polaronic phase is 
now very pronounced, being a consequence of (i). It is 
nonetheless interesting to observe that the critical density 
for the stabilization of the polaronic phase (rig ~ 0.004) 
and the typical stabilization temperature, Tm, are almost 
exactly the same as the ones encountered in Sec. III. 

The nature of the transitions as the temperature is 
lowered is as follows (see Fig. 16(b) for reference). Be- 
low TT-e ~ 0.001 the PM-Pol transition is continuous; 
np{T) varies continuously from to the saturation value 
Up = Tig, and no FM is stabilized except at very low 
temperatures, below the Pol phase. For Tig > 0.001, 
FM is stabilized with a continuous PM~FM transition; 
FM persists only until Tm{ne) is reached, at which point 
the polarons set in; the FM-Pol transition is discontin- 
uous because the magnetization drops to zero and Up 
jumps from to quasi-saturation upon crossing Tm- This 
can be seen clearly in the bottom frames of Fig. 16(b): 
above rif. > 0.001 the curve r,„ jumps discontinuously at 
Tmirie) to meet n'p. At lower temperatures, when the Tc 
line (red/circles) is crossed, the magnetization jumps to 
the value found in an homogeneous FM phase, concur- 
rently with a discontinuous drop of Hp from rig to zero. 
Notice that the Tc line is barely changed by the con- 
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FIG. 16: Polaronic stability relative to homogeneous 
PM/FM phases when a finite band filling is considered, (a) 
Phase diagram, (b) Each frame shows the equilibrium values 
of rip and M (order parameters) for total electronic concen- 
trations of ne = 0.0006, 0.0008, 0.0014, 0.0016. The blue (tri- 
angles) curve shows the equilibrium rip that is obtained when 
only paramagnetism is allowed (i.e. constraining M = 0). 



sideration of a finite band filling. This is related to the 
fact that, when Tc is reached, the polaron density has 
long since saturated at rig, thus emptying the band from 
carriers. 
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APPENDIX B: FINITE SIZE CORRECTIONS TO 
THE ELECTRONIC DOS 



with Hi G Z, which relates to ^{E) defined in (B2) by 



To estimate the 1/L corrections to the ground state 
energy of an electron gas consider free electrons inside 
a box of dimensions , Ly and . The electronic 
spectrum is 



2m 



with A; = TT ( Y-,^,Y- 

L/x J^y J^z 



The integrated DOS will clearly be 



n{E) 



E 

rij. .rij. >1 



{n, > 1) , 
(Bl) 

(B2) 



which corresponds, geometrically, to the set of integers 
{nx^ny^nz) bounded by the ellipsoid 



2 — T2 ^ T2 



TT 



LI 



(B3) 



In the thermodynamic limit, L^^y^z —>- oo and the usual 
procedure consists in replacing the discrete number of 
states satisfying (B3) by the volume of the ellipsoid in 
the first octant, divided by the elementary phase space 
volume. It is obvious that, by doing this, one is either 
neglecting or overcounting some of the points in phase 
space that rightfully satisfy (B3). This happens mainly 
at the boundaries, and is all right when L^.y^z ^ oo be- 
cause the errors arc of the order of 1/L or 1/i^ (Fig- 17). 
However, when L is finite, such corrections arc clearly 
relevant. A particularly important one arises from the 
fact that, when we calculate the volume of the ellipsoid, 
the points lying at the coordinate axes are automatically 
included, whereas from (Bl) they should not be. So let 



n'{E) = 



E 

TLj; .n^ ^riz 



9 (k^ - k{nx,ny,nzY) 



(B4) 
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(B5) 



The above result reflects the fact that, when calculating 
the continuum volume enclosed by the ellipsoid, one is 
adding an extra portion of phase space near the coor- 
dinate axes that should not be included, as in Fig. 17. 
Assuming that k^/tt^ is still reasonably greater than 
Lx ^,L~^ or L^^ , these terms correspond to 



Q! {E) ~ Vol. ellipsoid with axes Lx, Ly, Lz ; 
8 (k^ — k{ni, nj, 0)^) ~ Area of ellipse of axes Li, Lj ; 

ni.rij 

^ 6 (k^ - k{ni, 0, 0)^) ~ Length of the axis Li . 

(B6) 



Therefore, we have that 



4 f f k\^ 2 

-TT - LxLyLz = ?>n{E) + TT[-] {LxLy + LyLz + LzLx)-2-{Lx + Ly + Lz) + l. 

6 \Tr / \7I'/ TT 



(B7) 



implying that the corrected phase space volume is actually 



hi t\j 1 

^(E) = -^-^LxLyLz - —{LxLy + LyLz + LzLx) + '^i^x + Ly + Lz) — - ■ 



r 



(B8) 



There is still the error associated with the under/over 
estimates of the volume near the surface of the ellipsoid. 
This gives an additional contribution of the order C'(k)^^. 
But since the precise numerical factors are impossible to 
extract analytically, we do not include this correction. 
Consequently, the expression above is meaningful only 



down to ©(k^). Given that the volume and surface area 
of the volume enclosing the electron gas are related to the 
i's by LxLyLz = V and LxLy + LyLz + LzLx = S/2, 
the above can be summarized as 
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FIG. 17: Discrete Phase Space. 



from which the corrected DOS foUows: 
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For instance, taking the leading correction for electrons 
inside a cubic box {Li = L), implies a correction to the 
ground state energy per electron reading 
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with 
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In the main text, we are interested in the corrections 
to the energy of an electron gas confined to finite sized 
spherical bubbles. We notice that for a cube of side L, 
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(B13) 



inscribed spiicre 



the surface to volume ratio, equals the same ratio for the 
inscribed sphere with R = L/2. Hence we take the result 
(Bll) and simply substitute L by 2i?, obtaining 



E 



E, 



oo 



15 / TT \l/3 1 

16 V6n/ R 



(B14) 



which provides an estimate of the leading corrections to 
the ground state (and T = 0) energy of the confined 
electron gas. 
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